﻿using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Nature.Numerics;
using Nature.Toolkit;

namespace Nature.Data
{
    public class AH2004IdealGasSpeciesMolecularTransportProperties : IIdealGasSpeciesMolecularTransportProperties
    {
        #region Nested Types

        /// <summary>
        /// A set of indexes 
        /// </summary>
        struct IndexSet
        {
            /// <summary>
            /// 
            /// </summary>
            public int EpsilonIndex;

            /// <summary>
            /// 
            /// </summary>
            public int SigmaIndex;

            /// <summary>
            /// 
            /// </summary>
            public int MuIndex;

            /// <summary>
            /// 
            /// </summary>
            public int ReducedDipoleMomentIndex;
        }

        #endregion

        #region Private Fields

        private GasKinParameters[] m_gasKinValues;
        private IndexSet[][] m_binaryIndexes;

        private IObjectPool<double[]> m_speciesArrayPool;
        private IObjectPool<double[]> m_speciesLJPArrayPool;
        private IObjectPool<double[]> m_binaryLJPArrayPool;

        double[] m_binaryEpsilonDistinctValues;
        double[] m_binarySigmaDistinctValues;
        double[] m_binaryMuDistinctValues;
        double[] m_binaryReducedDipoleMomentDistinctValues;
        double[][] m_binaryDiffusivityConstantParts;
        double[][] m_binaryThermalDiffusivityConstantParts;

        double[] m_speciesEpsilonDistinctValues;
        double[] m_speciesReducedDipoleMoment;
        double[] m_speciesViscosityConstants;
        double[] m_speciesSelfDiffusivityConstantParts;
        int[] m_speciesEpsilonDistinctValueIndexes;

        #endregion

        #region Constructors

        public AH2004IdealGasSpeciesMolecularTransportProperties(AH2004ModelObject model)
        {
            this.Model = model;
            int nsp = model.Nsp;
            m_speciesArrayPool = ObjectPool<double[]>.CreateObjectPool(() => new double[nsp]);

            var w = model.SpeciesMolarMasses;

            this.Nsp = model.Nsp;
            m_gasKinValues = model.GasKinParameters.Data;

            IndexBuilder<double> epsilonIndexBuilder = new IndexBuilder<double>();
            IndexBuilder<double> sigmaIndexBuilder = new IndexBuilder<double>();
            IndexBuilder<double> muIndexBuilder = new IndexBuilder<double>();
            IndexBuilder<double> reducedMomentIndexBuilder = new IndexBuilder<double>();
            m_binaryIndexes = new IndexSet[Nsp][];
            m_binaryDiffusivityConstantParts = new double[Nsp][];
            m_binaryThermalDiffusivityConstantParts = new double[Nsp][];
            for (int isp = 0; isp < Nsp; ++isp)
            {
                m_binaryDiffusivityConstantParts[isp] = new double[isp + 1];
                m_binaryThermalDiffusivityConstantParts[isp] = new double[isp];
                IndexSet[] row = new IndexSet[isp + 1];
                m_binaryIndexes[isp] = row;
                GasKinParameters i = m_gasKinValues[isp];
                for (int jsp = 0; jsp <= isp; ++jsp)
                {
                    GasKinParameters j = m_gasKinValues[jsp];
                    GasKinParameters ij = new GasKinParameters(i, j);
                    double mass = (w[isp] * w[jsp] / (w[isp] + w[jsp])) / Constants.Avogadro; // molecular mass [kg]
                    double sig = ij.Sigma * Constants.Angstrom; // collision diameter [m]
                    double eps = ij.Epsilon * Constants.Boltzmann; // energy well depth [Joule]
                    double mu = ij.Mu * Constants.Debye; //	( Newton * meter ** 4) ** 0.5
                    m_binaryDiffusivityConstantParts[isp][jsp] =
                        3.0 / 16.0 * Math.Sqrt(2.0 * Math.PI * Math.Pow(Constants.Boltzmann, 3.0) / mass)
                        / Math.PI / (sig * sig) / Constants.Atmosphere;
                    double delta = 0.5 * (mu * mu) / eps / Math.Pow(sig, 3.0);
                    row[jsp] = new IndexSet
                    {
                        EpsilonIndex = epsilonIndexBuilder.Register(ij.Epsilon),
                        SigmaIndex = sigmaIndexBuilder.Register(ij.Sigma),
                        MuIndex = muIndexBuilder.Register(ij.Mu),
                        ReducedDipoleMomentIndex = reducedMomentIndexBuilder.Register(delta)
                    };
                    if (isp != jsp)
                    {
                        m_binaryThermalDiffusivityConstantParts[isp][jsp] = (w[isp] - w[jsp]) / (w[isp] + w[jsp]);
                    }
                }
            }
            m_binaryEpsilonDistinctValues = epsilonIndexBuilder.GetValues();
            m_binarySigmaDistinctValues = sigmaIndexBuilder.GetValues();
            m_binaryMuDistinctValues = muIndexBuilder.GetValues();
            m_binaryReducedDipoleMomentDistinctValues = reducedMomentIndexBuilder.GetValues();
            int binEpsLength = m_binaryEpsilonDistinctValues.Length;
            m_binaryLJPArrayPool = ObjectPool<double[]>.CreateObjectPool(() => new double[binEpsLength]);


            epsilonIndexBuilder = new IndexBuilder<double>();
            m_speciesEpsilonDistinctValueIndexes = new int[Nsp];
            m_speciesViscosityConstants = new double[Nsp];
            m_speciesReducedDipoleMoment = new double[Nsp];
            m_speciesSelfDiffusivityConstantParts = new double[Nsp];
            for (int isp = 0; isp < Nsp; ++isp)
            {
                GasKinParameters i = m_gasKinValues[isp];
                m_speciesEpsilonDistinctValueIndexes[isp] = epsilonIndexBuilder.Register(i.Epsilon);
                double mass = w[isp] / Constants.Avogadro; // molecular mass [kg]
                double sig = i.Sigma * Constants.Angstrom; // collision diameter [m]
                double eps = i.Epsilon * Constants.Boltzmann; // energy well depth [Joule]
                double mu = i.Mu * Constants.Debye; //	( Newton * meter ** 4) ** 0.5
                m_speciesViscosityConstants[isp] =
                    5.0 / 16.0 * Math.Sqrt(Math.PI * mass * Constants.Boltzmann) / Math.PI / sig / sig;
                m_speciesSelfDiffusivityConstantParts[isp] =
                    3.0 / 16.0 * Math.Sqrt(4.0 * Math.PI * Math.Pow(Constants.Boltzmann, 3.0) / mass)
                        / Math.PI / (sig * sig) / Constants.Atmosphere;
                m_speciesReducedDipoleMoment[isp] = 0.5 * (mu * mu) / eps / sig / sig / sig;
            }
            m_speciesEpsilonDistinctValues = epsilonIndexBuilder.GetValues();
            int ljpLength = m_speciesEpsilonDistinctValues.Length;
            m_speciesLJPArrayPool = ObjectPool<double[]>.CreateObjectPool(() => new double[ljpLength]);
        }

        #endregion

        #region Properties

        public ModelObject Model { get; private set; }

        protected int BinaryDistinctEpsCount { get { return m_binaryEpsilonDistinctValues.Length; } }

        #endregion

        #region Protected Virtual Methods

        protected virtual void CalcDiffusionDistinctLjp11(double temperature, double[] ljp11)
        {
            for (int i = 0; i < ljp11.Length; ++i)
            {
                ljp11[i] = Constants.LJP11(temperature / m_binaryEpsilonDistinctValues[i]);
            }
        }

        protected virtual void CalcThermDiffusionDistinctTauFactors(double temperature, double[] tauFactors)
        {
            using (var ljp11OpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            using (var ljp22OpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            using (var ljp12OpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            using (var ljp13OpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            using (var tauOpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            using (var reducedTemperatureOpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            {
                double[] ljp11 = ljp11OpjectPoolItem.Unwrap();
                double[] ljp22 = ljp22OpjectPoolItem.Unwrap();
                double[] ljp12 = ljp12OpjectPoolItem.Unwrap();
                double[] ljp13 = ljp13OpjectPoolItem.Unwrap();
                double[] reducedTemperature = reducedTemperatureOpjectPoolItem.Unwrap();
                Debug.Assert(reducedTemperature.Length == tauFactors.Length);
                for (int i = 0; i < reducedTemperature.Length; ++i)
                {
                    reducedTemperature[i] = temperature / m_binaryEpsilonDistinctValues[i];
                    ljp11[i] = Constants.LJP11(reducedTemperature[i]);
                    ljp22[i] = Constants.LJP22(reducedTemperature[i]);
                    ljp12[i] = Constants.LJP12(reducedTemperature[i]);
                    ljp13[i] = Constants.LJP13(reducedTemperature[i]);

                    double aij = ljp22[i] / ljp11[i];
                    double bij = (5.0 * ljp12[i] - 4.0 * ljp13[i]) / ljp11[i];
                    double cij = ljp12[i] / ljp11[i];

                    tauFactors[i] = (15.0 / 2.0) *
                        (2.0 * aij + 5.0) * (6.0 * cij - 5.0) /
                        (aij * (16.0 * aij - 12.0 * bij + 55.0));
                }
            }
        }

        #endregion

        #region IIdealGasSpeciesMolecularTransportProperties Members

        public int Nsp { get; private set; }

        public virtual void CalcSpeciesIdealGasViscosities(double temperature, double[] viscosities)
        {
            double tSqrt = Math.Sqrt(temperature);
            using (var ljp22ArrayObjectPoolItem = m_speciesLJPArrayPool.Acquire())
            using (var tsArrayObjectPoolItem = m_speciesLJPArrayPool.Acquire())
            {
                double[] ljp22 = ljp22ArrayObjectPoolItem.Unwrap();
                double[] ts = tsArrayObjectPoolItem.Unwrap();
                for (int i = 0; i < ljp22.Length; ++i)
                {
                    ts[i] = temperature / m_speciesEpsilonDistinctValues[i];
                    ljp22[i] = Constants.LJP22(ts[i]);
                }

                for (int isp = 0; isp < Nsp; ++isp)
                {
                    int ix = m_speciesEpsilonDistinctValueIndexes[isp];
                    double om22 = ljp22[ix];
                    double delta = m_speciesReducedDipoleMoment[isp];
                    om22 += 0.2 * delta * delta / ts[ix];
                    viscosities[isp] = m_speciesViscosityConstants[isp] * tSqrt / om22;
                }
            }
        }

        public virtual void CalcSpeciesIdealGasThermalConductivities(double temperature, double[] conductivities)
        {
            var w = Model.SpeciesMolarMasses;
            var thermo = Model.IdealGasSpeciesThermoPotentials;

            using (var viscositiesObjectPoolItem = m_speciesArrayPool.Acquire())
            using (var selfDiffusivitiesObjectPoolItem = m_speciesArrayPool.Acquire())
            using (var speciesCvObjectPoolItem = m_speciesArrayPool.Acquire())
            {
                double[] viscosities = viscositiesObjectPoolItem.Unwrap();
                double[] selfDiffusivities = selfDiffusivitiesObjectPoolItem.Unwrap();
                double[] speciesCv = speciesCvObjectPoolItem.Unwrap();

                this.CalcSpeciesIdealGasViscosities(temperature, viscosities);
                this.CalcSpeciesIdealGasSelfDiffusivities(temperature, Constants.Atmosphere, selfDiffusivities);
                thermo.Cv(temperature, speciesCv);


                double Ctrans, Crot, Cvib;
                double Ftrans, Frot, Fvib;
                double ro0, A, B, Zrot;
                double cond;

                for (int isp = 0; isp < Nsp; ++isp)
                {
                    var gasKin = m_gasKinValues[isp];
                    double visc = viscosities[isp];
                    double Dkk = selfDiffusivities[isp];
                    double cv = speciesCv[isp];
                    //		

                    ro0 = w[isp] * (Constants.Atmosphere / (Constants.Rgas * temperature));
                    //
                    Zrot = gasKin.Zrot *
                        gasKin.RotRelCollFactor(298.0) /
                        gasKin.RotRelCollFactor(temperature);

                    // 
                    switch (gasKin.GeometryIndex)
                    {
                        case (0)://	single atom
                            {
                                Ctrans = 1.5;
                                Crot = 0.0;
                                Cvib = 0.0;
                            }
                            break;
                        case (1)://	linear molecule
                            {
                                Ctrans = 1.5;
                                Crot = 1.0;
                                Cvib = Math.Max(0.0, cv - 2.5);
                            }
                            break;
                        case (2)://	nonlinear molecule
                            {
                                Ctrans = 1.5;
                                Crot = 1.5;
                                Cvib = Math.Max(0.0, cv - 3.0);
                            }
                            break;
                        default:
                            throw new InvalidOperationException();
                    }//switch( part.geometry )
                    //
                    //
                    double factor = ro0 * Dkk / visc;
                    A = 2.5 - factor;
                    B = Zrot + 2.0 / Math.PI * (5.0 / 3.0 * Crot + factor);
                    //
                    Ftrans = 2.5 * (1.0 - 2.0 / Math.PI * Crot / Ctrans * A / B);
                    //
                    Frot = factor * (1.0 + 2.0 / Math.PI * A / B);
                    //
                    Fvib = factor;
                    //
                    cond = visc / w[isp];
                    cond *= (Ftrans * Ctrans + Frot * Crot + Fvib * Cvib);
                    cond *= Constants.Rgas;
                    conductivities[isp] = cond;
                }
            }
        }

        public virtual void CalcSpeciesIdealGasSelfDiffusivities(double temperature, double pressure, double[] selfDiffusivities)
        {
            double reducedPressure = pressure / Constants.Atmosphere;
            double t3Sqrt = Math.Pow(temperature, 3.0 / 2.0);
            using (var ljpArrayObjectPoolItem = m_speciesLJPArrayPool.Acquire())
            using (var tsArrayObjectPoolItem = m_speciesLJPArrayPool.Acquire())
            {
                double[] ljp11 = ljpArrayObjectPoolItem.Unwrap();
                double[] ts = tsArrayObjectPoolItem.Unwrap();
                for (int i = 0; i < ljp11.Length; ++i)
                {
                    ts[i] = temperature / m_speciesEpsilonDistinctValues[i];
                    ljp11[i] = Constants.LJP11(ts[i]);
                }

                for (int isp = 0; isp < Nsp; ++isp)
                {
                    int ix = m_speciesEpsilonDistinctValueIndexes[isp];
                    double om11 = ljp11[ix];
                    double delta = m_speciesReducedDipoleMoment[isp];
                    om11 += 0.19 * delta * delta / ts[ix];
                    selfDiffusivities[isp] = m_speciesSelfDiffusivityConstantParts[isp] * t3Sqrt / om11 / reducedPressure;
                }
            }
        }

        public void CalcSpeciesIdealGasDiffusivities(double temperature, double pressure, double[][] diffusivities)
        {
            double reducedPressure = pressure / Constants.Atmosphere;
            double t3Sqrt = Math.Pow(temperature, 3.0 / 2.0);
            using (var ljpOpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            using (var reducedTemperatureOpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            {
                double[] ljp11 = ljpOpjectPoolItem.Unwrap();
                double[] reducedTemperature = reducedTemperatureOpjectPoolItem.Unwrap();
                this.CalcDiffusionDistinctLjp11(temperature, ljp11);
                for (int i = 0; i < reducedTemperature.Length; ++i)
                {
                    reducedTemperature[i] = temperature / m_binaryEpsilonDistinctValues[i];
                }
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    double[] disp = diffusivities[isp];
                    double[] constParts = m_binaryDiffusivityConstantParts[isp];
                    for (int jsp = 0; jsp <= isp; ++jsp)
                    {
                        IndexSet ix = m_binaryIndexes[isp][jsp];
                        double tstar = reducedTemperature[ix.EpsilonIndex];
                        double om11 = ljp11[ix.EpsilonIndex];
                        double delta = m_binaryReducedDipoleMomentDistinctValues[ix.ReducedDipoleMomentIndex];
                        om11 += 0.19 * delta * delta / tstar;
                        disp[jsp] = constParts[jsp] * t3Sqrt / om11 / reducedPressure;
                        Debug.Assert(disp[jsp] > 0.0d && disp[jsp] < 1.0d,
                            string.Format("[{0}][{1}]", isp, jsp));
                    }
                }
            }
        }

        public void CalcSpeciesIdealGasThermoDiffusivities(double temperature, double[][] thermDiffusivities)
        {
            using (var tauOpjectPoolItem = m_binaryLJPArrayPool.Acquire())
            {
                double[] tau = tauOpjectPoolItem.Unwrap();
                this.CalcThermDiffusionDistinctTauFactors(temperature, tau);
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    double[] disp = thermDiffusivities[isp];
                    double[] constParts = m_binaryThermalDiffusivityConstantParts[isp];
                    IndexSet[] ix = m_binaryIndexes[isp];
                    for (int jsp = 0; jsp < isp; ++jsp)
                    {
                        int index = ix[jsp].EpsilonIndex;
                        disp[jsp] = tau[index] * constParts[jsp];
                    }
                }
            }
        }

        #endregion
    }

    public class AH2004IdealGasSpeciesMolecularTransportPropertiesTabulated : AH2004IdealGasSpeciesMolecularTransportProperties
    {
        public readonly double Tmin;

        public readonly double Tmax;

        public readonly double Delta;

        #region Private Fields
        private Tabulated m_viscosities;
        private Tabulated m_conductivities;
        private Tabulated m_selfDiffusivities;
        private Tabulated m_thermDiffTau;
        private Tabulated m_binaryLjp11;
        #endregion

        public AH2004IdealGasSpeciesMolecularTransportPropertiesTabulated(AH2004ModelObject model, double tmin, double tmax, double delta)
            : base(model)
        {
            Tmin = tmin;
            Tmax = tmax;
            Delta = delta;
        }

        public override void CalcSpeciesIdealGasViscosities(double temperature, double[] viscosities)
        {
            if (m_viscosities == null)
            {
                VectorFunctionDelegate func = (_visc, _temp) => { base.CalcSpeciesIdealGasViscosities(_temp, _visc); };
                m_viscosities = new Tabulated(func, Nsp, Tmin, Tmax, Delta);
            }
            m_viscosities.Approximate(viscosities, temperature);
        }

        public override void CalcSpeciesIdealGasThermalConductivities(double temperature, double[] conductivities)
        {
            if (m_conductivities == null)
            {
                VectorFunctionDelegate func = (_cond, _temp) => { base.CalcSpeciesIdealGasThermalConductivities(_temp, _cond); };
                m_conductivities = new Tabulated(func, Nsp, Tmin, Tmax, Delta);
            }
            m_conductivities.Approximate(conductivities, temperature);
        }

        public override void CalcSpeciesIdealGasSelfDiffusivities(double temperature, double pressure, double[] selfDiffusivities)
        {
            if (m_selfDiffusivities == null)
            {
                VectorFunctionDelegate func = (_diff, _temp) => { base.CalcSpeciesIdealGasSelfDiffusivities(_temp, Constants.Atmosphere, _diff); };
                m_selfDiffusivities = new Tabulated(func, Nsp, Tmin, Tmax, Delta);
            }
            m_selfDiffusivities.Approximate(selfDiffusivities, temperature);
            if (pressure != Constants.Atmosphere)
            {
                double reducedPressure = pressure / Constants.Atmosphere;
                for (int isp = 0; isp < Nsp; ++isp) { selfDiffusivities[isp] /= reducedPressure; }
            }
        }

        protected override void CalcThermDiffusionDistinctTauFactors(double temperature, double[] tauFactors)
        {
            if (m_thermDiffTau == null)
            {
                VectorFunctionDelegate func = (_tau, _temp) => { base.CalcThermDiffusionDistinctTauFactors(_temp, _tau); };
                m_thermDiffTau = new Tabulated(func, BinaryDistinctEpsCount, Tmin, Tmax, Delta);
            }
            m_thermDiffTau.Approximate(tauFactors, temperature);
        }

        protected override void CalcDiffusionDistinctLjp11(double temperature, double[] ljp11)
        {
            if (m_binaryLjp11 == null)
            {
                VectorFunctionDelegate func = (_ljp11, _temp) => { base.CalcDiffusionDistinctLjp11(_temp, _ljp11); };
                m_binaryLjp11 = new Tabulated(func, BinaryDistinctEpsCount, Tmin, Tmax, Delta);
            }
            m_binaryLjp11.Approximate(ljp11, temperature);
        }
    }
}
